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ABSTRACT 

Aims. We extend the gr-band time coverage of the gravitationally lensed double quasar Q0957+561. New gr light curves permit us 
to detect significant intrinsic fluctuations, to determine new time delays, and thus to gain perspective on the mechanism of intrinsic 
variability in Q0957-I-561. 

Methods. We use new optical frames of Q0957+561 in the g and r passbands from January 2005 to July 2007. These frames are part of 
an ongoing long-term monitoring with the Liverpool robotic telescope. We also introduce two photometric pipelines that are applied 
to the new gr frames of Q0957+561. The transformation pipeline incorporates zero-point, colour, and inhomogeneity corrections to 
the instrumental magnitudes, so final photometry to the 1-2% level is achieved for both quasar components. The two-colour final 
records are then used to measure time delays. 

Results. The gr light curves of Q0957+561 show several prominent events and gradients, and some of them (in the g band) lead 
to a time delay between components AtgA = 417 ± 2 d (Icr). We do not find evidence of extrinsic variability in the light curves of 
Q0957-I-561. We also explore the possibility of a delay between a large event in the g band and the corresponding event in the r band. 
The gr cross-correlation reveals a time lag At,.g = 4.0 ± 2.0 d (Icr; the g-band event is leading) that confirms a previous claim of the 
existence of a delay between the g and r band in this lensed quasar 

Conclusions. The time delays (between quasar components and between optical bands) from the new records and previous ones in 
similar bands indicate that most observed variations in Q0957+561 (amplitudes of ~ 100 mmag and timescales of ~ 100 d) are very 
probably due to reverberation within the gas disc around the supermassive black hole. 

Key words, techniques: photometric - gravitational lensing - black hole physics - quasars: individual: Q0957+561 



1. Introduction 

Studies of optical continuum variability in gravitationall y lensed quasars (GLQs) have a rnain advantage: one is usually able to 
disentangle intrinsic from extrinsi c signal in GLQ s (e.g.. lKundic et al.1ll997l:|Paraficz et al. Il2006t iGoicoechea et aLll2008l Paper 
I). Following the original idea by iRefsdall ( 1 19641) . intrinsic variations in brightness records of GLQs have mainly been used to 
estimate globa l time delays between components, and to discuss the structure of galaxy mass halos and the expansion rate of the 
Universe (e.g.. lKochanek et al. 1 12004, and references therein). Less effort has been devoted to investigating the nature of intrinsic 
fluctuations, which are generated by mechanisms of variability in lensed quasars. This can be done by measuring time delays 
between components and between optical bands, using prominent events in segments of long-term light curves. Time delays between 
two given components of a GLQ (determined from different pairs of twin features) arise from gravitational lensing of flares in the 
variable source. While the flaring of a well-defined emission region (e.g., a ring of the accretion disc) produces a set of similar 
delays, the existence of flares in some widely separated zones can lead to important time delay differences (lYonehara I[l999 ). For 
either of the two components, time delays between optical bands (or interband time delays) refer to time lags arising from physical 

phenomena within the quasar. 

The gravitationally lensed double quasar Q0957H-561 at z = 1.41 (IWalsh et al. Ill979h has been monitored photometrically in 
different optical bands with different telescopes. To derive a global time delay between quasar components, some previous studies 
used large data sets incorporating all kinds of fluctuations, i.e., noisy or poorly sampled features as well as noticeable gradients and 
events on several timescales. These large data sets are based on frames that were taken in the 1980s and 90s, and they lead to a global 
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delay of about 423-425 d dOscoz et al. 11200 lllOvaldsen et al. Il2003ah . The Apache Point Observatory (APO) experi ment permitted 
investigators to follow-up the variability in the g and r bands during the 1995 and 1996 seasons, i.e., covering 1 .5 years ("Kundi c et al. I 
[1997). This monitoring programme with the APO 3.5 m telescope produced accurate light curves of both components Q0957-I-561 A 
and Q0957+561B, which show sharp intrinsic features with high signal-to-n oise ratio (S /N > 3) . From the APO main twin events 
in the g band (a prominent event in A and the replica event in B; S /N ~ 6.5). lKundic et al. 1 (ll997l) also reported a gravitational lens 
time delay of 417 '^^ d (95% confidence interval). Complementary to this result. Collier I (1200 lb found that the r-band main twin 
events lag with respect to the ones in the ^-band by 3.4 d (68% confidence interval), and this interband delay was interpreted as 
clear evidence for reprocessing in the accretion disc of the quasar. 

iGoicoechea I ( l2002b reanalysed the APO g-hand data set to obtain two different gravitational lens time delays of 417.0 + 0.6 d 
(68% confidence interval) and 432.0 + 1.9 d (68% confidence interval) depending on the features taken as reference. The longest 
delay corresponds to the APO secondary twin events {S /N ~ 3) and it clearly disagrees with the 417-day value. The APO main 
and secondary twin events in the g band are associated with a main flare and a secondary flare in the variable source, respectively. 
From the time delay difference of 15 + 2 days (68% confidence interval) and using standard cosmological parameters (Q - 0.3, A - 
0.7), one can also determine a minimum size fo r the variable source (minimum distance between both flares) of 300 pc dYoneharal 
1 1 9991: IGoicoechea II2002I: [Yonehara et al. 1 12003). Several physical sources are consistent with this spatial constraint (from multiple 
gravitational lens delays) and the interband delay for the main twin events, but either a nuclea r accretion disc pl us a circumnuclear 
stellar region or a nuclear accretion disc plus an optical jet are the most probable ones (e.g., lHutchingsll2003l) . We note that the 
424-d global time delay between components dOscoz et al. 1 1200 lUOvaldsen et al. Il2003ah coincides with the average of both APO 
gravitational lens time delay s. This suggests the presence of two or more gravitational lens delays in the current large data sets (for 
alternative explanations, see lSchild |[2005l: iHirv et al. Il2007h . 

In the present study, we substantially extend the gr-hand time coverage of Q0957+56 1 . The key idea is to use new gr light curves 
to detect prominent intrinsic events similar to the APO ones. These new features should allow us to determine new time delays and 
to improve our understanding of the mechanism causing the intrinsic variability. The paper is organised as follows: in Section 2, we 
present new data of Q0957-t-561 based on recent observations with the Liverpool 2 m telescope (LT) in the g and r bands, spanning 
2.5 years. We describe the observations, the pre-processing, and the photometric procedure for determining calibrated and corrected 
magnitudes of field stars and quasar components. This last reduction procedure consists of two new pipelines specially designed for 
the LT. In Section 3 we study the time delays between the two components of Q0957H-561 as well as the possible delays between 
the g and r band in the new data set. In Section 4 we summarize our results. From the APO and LT delays of Q0957+561, we also 
discuss the origin of the observed variations with an amplitude of ~ 100 mmag and lasting ~ 100 d. 



2. Data acquisition and reduction 

2.1. Observations and pre-processing 

Liverpool Quasar Lens Monitoring (LQLM) I is the first phase o f an optical follow- up of lensed quasars, undertaken using the 
RATCam optical CCD camera on the Liverpool robotic telescope jSteele et al. 112004b between January 2005 and July 2007. The 
first scientific output of LQLM I was reported in Paper I, and we concentrate here on the observations of Q0957-)-561 in the g and 
r filters. The field of view and the pixel scale (binning 2x2) were ~ 4'6 x 4'6 and 0'.'278, respectively. The exposure times were 
100-200 s (g band) and 120 s (r band). We obtained 286 frames in the g band and 264 frames in the r band. The LT observed for a 
total science time of ~ 22.6 h during the 2.5-year programme of Q0957-I-561. 

A pre-processing pipeline is applied to all RATCam framefl This performs three basic instrumental reductions: bias subtraction, 
trimming of the overscan regions, and flat fielding. We also apply a bad-pixel mask (made available by the Angstrom project; 
[Kerins et al. 2006), and correct bad pixels on the CCD. The next step is the pre-selection of frames, based on individual inspection, 
to assure that exposures verify some elemental conditions (e.g., that the telescope pointing was accurate enough so that the lens 
system was included in the field of view, that there is no strongly degraded signal, etc), and that seeing {FWHM) and sky level 
(background) values do not exceed reasonable bounds. We only consider frames with FWHM < 3" due to the separation between 
the two quasar components (A and B) of ~ 6". The pre-selected database contains 199 frames in the g band and 210 frames in the r 
band. This means that ~ 75% of the original LT frames were initiaUy useful. 



2.2. Instrumental photometry 

In a first step, we take a reference frame, i.e., a high-quality frame with small FWHM and large signal-to-noise ratio (SNR). We 
then measure the positions (with respect to the left bottom corner of the reference frame )of seven reference stars and b oth quasar 
images. We select the 7 brightest stars in the Sloan Digital Sky Survey (SDSS) catalog^ (e.g., Adelman-M cCarthv et al. 2007]). 
These stars, having ^(SDSS) and r(SDSS) magnitudes below 18 and 17, respectively, were labeled as X, G, F, H, D, E, and R stars 



' See the Web site http://telescope.livjm.ac.uk/Info/TelInst/Inst/RATCam/index.phpl 

^ See the DR6 Catalogue Archive Server site http://cas.sdss.org/astrodr6/en/ Funding for the SDSS has been provided by the Alfred R Sloan 
Foundation, the Participating Institutions, the NASA, the NSF, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck 
Society. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions 
are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los 
Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico 
State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. 
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Fig. 1. Observations of Q0957+561 with the Liverpool robotic telescope in the g band. We display system subframes (left panels), 
model subframes (middle panels), and residual subframes (right panels) for five frames taken during the 2.5-year monitoring period 
(see main text). 

in Figure 1 and Table 1 of lOvaldsen et al. I ( l2003al) . Several Image Reduction and Analysis Facility (IRAFj^ tasks are also used to 
identify the available reference objects (in general, less than 7 stars) and the quasar components in the rest of the frames. 

In a second step, our photometric pipeline performs aperture photometry of bright field stars and quasar images. This IRAF 
procedure is used to estimate initial instrumental fluxes (sources and their associated backgrounds) and to improve the initial source 
positions on each frame. The pipeline also cuts the original frames in order to produce square subframes with 64 pixels per side: 
the system subframe (around the centre of the lens system) and subframes of stars (around the bright stars), and makes a PSF 
subframe containing the clean 2D profile of the H star (removing the local background). This last empirical PSF is required when 
performing PSF fitting. The point-like sources (quasar components and stars) are modelled by means of the empirical PSF, whereas 
the extended source (lensing elliptical galaxy) is modelled by a de Vaucouleurs profile convolved with the empirical PSF. After 
obtaining all subframes for a given frame, PSF photometry on the stellar and system (crowed field) subframes is performed with 
IMFITFITS software ( McLeod et al. 1998). The pipeline is written in the Python programming languag^ and incorporates the 
capabilities of IRAF (through the PyRAF interface) and IMFITFITS, as well as additional numerical and graphical tools. 

To determine accurate quasar fluxes, one needs to use a set of constraints. F or Q0957H-561 , the most relevant constraints were 
obtained from Hubble Space Telescope (HST) frames dBernstein e t al. 1997; Ke eton et al. 11199 8; Kochanek et al. 2008): positions 
of the B component and the lensing galaxy relative to the A component, and the optical structure of the galaxy, i.e., eff'ective radius, 
ellipticity, and position angle (a de Vaucouleurs profile was fitted to HST images). Due to the relatively low brightness of the lensing 
galaxy in the frames and the proximity of the B component to the galaxy, we determine the galaxy-to-H star ratio (GAL/H) in the 
gr bands from the best LT frames, in terms of FWHM and SNR values. The H star is relatively bright and it is present in all frames. 
We then apply the pipeline to all frames (whatever their qualities) in each optical filter, by setting the galaxy fluxes to those derived 
from the GAL/H ratio and H star fluxes, and allowing the remaining free parameters to vary. 

The photometry pipeline output includes the system subframes, their model subframes (best fits) and the associated residual 
subframes (system subframes after subtracting model subframes). In Fig. [1] we show system subframes (left panels), model sub- 
frames (middle panels), and residual subframes (right panels) corresponding to five frames in the g band. From top to bottom: 
March 16, 2005 (FWHM = 2'.'03, SNR = 215, x' = 1-44), November 9, 2005 (FWHM = 1'.'75, SNR = 313, = 1-32), April 

^ IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in 
Astronomy (AURA) under cooperative agreement with the National Science Foundation. 
See the Web site http://www.python.org/ , 
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Fig. 2. Colour coefficient in the g band. The values are distributed around the central discontinuous line (average coefficient), and 
most of them are placed between the top and bottom discontinuous lines (filled circles). Only seven extreme values (triangles and 
open circles) exceed these limits. 



26, 2006 (FWHM = 1'.'76, SNR = 371, = 2.15), October 21, 2006 (FWHM = 1'.'48, SNR = 321, = 1-30), and May 28, 
2007 (FWHM - r.'63, SNR = 120, x^ - 1-43), where the SNR values ai-e infeiTed from the A images (having fluxes similar to 
those of B images) and the;i'^ values quantitatively describe the quality of the fits (i.e., these represent the standard reduced for 
the best fits). All subframes in Fig. [T]have been expanded by a factor of 2. The visual comparison between left and middle panels 
as well as the inspection of patterns of residual brightness (right panels) indicate that the photometric method works well. The 
pipeline also produces a basic data release file containing values of all relevant instrumental fluxes (stars and quasar images) and 
relative instrumental magnitudes of both quasar components, e .g., ^* - g*^ and e* - g^ in the g band. To check the reliability of our 
PSF fitting procedure, we applied a deconvolution technique jKoptelova et al. 1 12005) to two sets of frames in October-December 
2005 igr bands). The re lative instrumental magn itudes from the deconvolution method agreed with the records from the PSF fitting 
technique (see Fig. 3 of lGoicoechea et al7ll2007h . 

2.3. Calibrated and corrected magnitudes 

We use a transformation pipeline (in the Python programming language) to obtain SDSS magnitudes from instrumental magnitudes 
that are corrected for systematic eff'ects. The whole calibration-correction process is outlined in Appendix A. Only frames with 
SNR > 100 over Q0957-I-561A are taken into account (see however Section 3). In the g band, this selection leads to 170 frames. 
In the r band, besides the SNR based selection, the surviving frames from the first season (January-June 2005) are also removed. 
Several of these ~ 10 r-band frames with SNR above 100 (first season) are characterized by an anomalous image formation. Thus, 
the high-quality data set in the r band includes 167 frames. 

The transformation pipeline fits the deviations between instrumental and standard g magnitudes of the 7 reference stars to the 
transformation model that incorporates a zero-point term (ag), a colour coefficient (J3g), and inhomogeneity coefficients (yg,„m). Eq. 
(A.l 1) shows the relationship between the observed magnitude deviation and the model to describe it. The zero-point term and the 
colour coefficient are allowed to vary over time because the atmospheric and instrumental conditions significantly evolve during 
the 2.5 years of monitoring. The last ingredient of the model is a linear-quadratic inhomogeneity term, which is r elated to the 2D 
position on the CCD and trie s to correct the possible inhomogeneous response over the camera area (e.g., .Manfroidetal.1120011: 
iMagnier & Cuillandre II2004I) . Each source occupies different positions on the CCD area during the robotic monitoring period, so 
this could complicate the collecting of accurate brightness records. 

With respect to the least squares fit, in Fig.|2]we plot the solution of /3g. The /3g values are distributed around an average colour 
coefficient (/3g} = -0.097 (central discontinuous line in Fig.|2]i, which is close to the typical coefficient (see comments in Appendix 
A). The scatter is cr(J3g) - 0.033, and the {J3g} + 2.5o-(J3g) limits also appear in Fig. |2] (top and bottom discontinuous lines). In 
relation to the average coefficient, there are seven extreme values representing changes from 100%, i.e., values around either -0.2 
or 0.0 (triangles and open circles in Fig. |2]i. The first two triangles and the first open circle correspond to the first night after a 
realuminisation of the mirror and other maintenance works, and very probably, the telescope was not performing optimally that 
night. The rest of the extreme values (around day 3770) are associated with dates close to periods of very bad weather The highest 
values of j3g are a consequence of atmospheric-instrumental perturbations during the hard winter in January-February 2006. From the 
best solutions of jg,„„,, the pipeline also produces the 2D inhomogeneity pattern, i.e., 2o<«+m<2 7 g,nmx"y'^ ■ This is depicted in Fig. [3] 
In the transformation procedure, we set the origin of coordinates at the centre of the 1024x1024 CCD, so it has an inhomogeneity 
correction equal to zero . In Fig.[3l we see an inhomogeneity amplitude of ~ 80 mmag, which is consistent with results from other 
optical telescopes (e.g.. iManfroid et al. II2001I : iMagnier & Cuillandre |[2004l) . It is evident that the inhomogeneity pattern in Fig.|3] 
plays a role in achieving 1 -2% photometric accuracy. 
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Fig. 3. Inhomogeneity map in the g band. The zero inhomogeneity level is described by means of a continuous line that crosses the 
centre of the 1024x1024 camera. Pixels inside this zero level contour have a positive inhomogeneity of a few mmag, whereas the 
rest of the pixels have a negative inhomogeneity that varies between a few mmag and tens of mmag. We explicitly show contours 
of -10, -20, -30, -40, -50, -60, and -70 mmag. 
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Fig. 4. Final magnitudes of Q0957+561A (top panel) and Q0957+561B (bottom panel) in the g band of the SDSS photometric 
system. These ^-SDSS light curves include noticeable fluctuations covering a 2.5-year monitoring period from January 2005 to 
June 2007. 



After the g-band fit, the pipeline computes the calibrated and corrected records of the reference stars and both quasar images 
from Eq. (A. 13). The 14-15th magnitude stars have a typical scatter of 5 mmag, the 15-16th magnitude stars are characterized by 
a slightly larger scatter of about 7 mmag, and the faintest ~ 18th magnitude star (R star) has a scatter of about 17 mmag. Although 
we have several nights with two or three exposures at different times, the standard intranight deviations of the stellar curves do not 
trace their scatters (see however Paper I). This is not surprising because the intranight variations exclusively correspond to several 
nights in the second season, which covers a small fraction of the total monitoring period. Thus, after some preliminary test using 
the stellar records, we find a non-biased estimator of uncertainties (typical errors): stellar scatters are well traced by the standard 
deviations between adjacent magnitudes that are separated by < 3 d. As this error estimator gives reasonable results, we apply it to 
the g-SDSS magnitudes of Q0957+561. 
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Fig. 5. Final magnitudes of Q0957+561A (top panel) and Q0957+561B (bottom panel) in the r band of the SDSS photometric 
system. The r-SDSS records from October 2005 to June 2007 (two whole seasons) incorporate different prominent features that are 
also seen in the g-SDSS curves (see Fig.|4]i, with the g-SDSS features having a larger amplitude. 



A final refinement (selection) is taken into account. Our last selection criterion is colour based: frames with extreme colour 
coefficients (see the triangles and open circles in Fig.|2]i are also removed from the data set. This leads to 163 surviving frames. 
We obtain uncertainties (see above) of about 16 mmag in both ~ 17th magnitude quasar components, i.e., photometry to the 1-2% 
is achieved for the lensed quasar These typical errors are in complete agreement with the stellar scatters, since they are clearly 
larger than 5-7 mmag (results for the brightest reference stars) and similar to the result for the faintest reference star R. For each 
component, we also group pairs or trios of magnitudes measured on the same night. The final light curves of Q0957+561A (top 
panel) and Q0957+561B (bottom panel) are shown in Fig.|4] These g-SDSS light curves include important gradients and prominent 
events, which resemble those reported bv lKundic et al. 1 (11997) using APO observations in the g band. The whole monitoring period 
consists of three observational seasons: January- June 2005 (first season), October 2005-June 2006 (second season), and October 
2006-June 2007 (third season). Besides the three observational seasons, there are two important gaps in the LT monitoring as a 
consequence of the annual occultation of the lens system. 

The whole procedure in the g band is repeated in the r band. All frames with extreme colour coefficients are not considered 
in building the final light curves, so we use a data set incorporating 142 frames. With respect to the quasar brightness records, 
we achieve ~ 1% photometry (errors of about 12 mmag). The final (grouped) magnitudes are presented in Fig.|5] where the top 
and bottom panels display the records of Q0957+561A and Q0957-I-561B, respectively. These r-SDSS light curved trace promi- 
nent fluctuations that are weaker tha n the corr esponding fluctuations in g-SDSS (see Fig. |4] and subsection 3.2). A similar result 
was claimed by th e APO team (Kun dic et al. I 1995, 19 97), and some evid ence for chromatic variability was also suggested by 
lUUan et alTI (l2003b (see also the BVRI variations in Ser ra-Ricart et al. Ill999l). The LT records in the red region of the optical spec- 
trum are less noisy than previous curves at red wavelengths (e.g.. lKundic et al. Ill997t ISerra-Ricart et al. 1 [l999l) . which is due to a 
combination of an absence of relatively short variations and strict selection procedures. 



3. Time delays of Q0957+561 

3.1. Delay between quasar components 

The ^-band light curve of A in the second season (October 2005-June 2006) shows significant fluctuations that are repeated in the 
;^- band light curve of B during the third season (October 2006-June 2007). Taking int o account the expected delay range of 415-435 
d (iKundic et al. lll997HOscoz et al. l200ltrGoicoecheal2'00llOvaldsen et al. Il2003ah . this result is fully consistent with the presence 
of intrinsic fluctuations in those records. We use the ^-SDSS magnitudes of A and B in the second and third seasons, respectively, 
to accurately measure the time delay (s) between both components of Q0957 -1-561. 



^ The gr records are available at |http://grupos.umcan.es/glendama/| 



Shalyapin et al.: New two-colour light curves of Q0957+561 



7 




MJD-50000 

Fig. 6. Comparison between the ^-band light curve of A in the second season (shifted by 420 d; see main text) and the g-hand light 
curve of B in the third season. The A record (filled circles) shows two different features separated by a gap of about 50 d: while the 
first feature contains an event AElg and the beginning of another consecutive event AE2g, the second feature describes the (noisy) 
decline in flux of AE2g. A vertical line is drawn to distinguish between the two events AElg and AE2g. Replica events BElg and 
BE2g are clearly seen in the B record (open circles). 



Table 1. Magnitude offset and time delay measurements in the g band. 



Brightness records 


Method 


Offset^ (mag) 


Delay_^ (d) 


A(season 2)-B(season 3) 
AEl-BEl 


x' 


-0.090 + 0.004 
-0.092 + 0.004 
-0.083 ± 0.006 


417 ±2 

416 ±5 

417 ± 2 
417 ± 2 



" Magnitude offset between the A and B components, where the sign "-" means that A is fainter (see Fig.|6](. From S' we do not measure the 
shift in magnitude, since 5^ is a technique based on autocorrelation and cross-correlation functions. All measurements are Itr intervals. 
* Delay of the replica variation in B with respect to the variation in A (the A component is leading). All measurements are Icr intervals. 



About one half of the frames with 80 < SNR < 100 coiTespond to the second season, and thus, they could help to trace the 
variability of A and to minimize uncertainties in time delay estimates. Their photometric outputs (magnitudes of A) are consistent 
with results from SNR > 100 frames at adjacent epochs, so we recover them and expand the g-band record of A in the second season. 
In Fig.|6] the A light curve, shifted by 420 d (filled circles), and the unchanged B light curve (open circles) are plotted. A reference 
value of 420 d is used to shift in time one component and to compare it with the other (see above and Introduction). The A record 
shows two different features separated by a gap of about 50 d (caused by atmospheric-instrumental problems in January-February 
2006; see subsection 2.3). The first feature in the A curve consists of an event AElg and the beginning of another consecutive event 
AE2g, whereas the second feature is a noisy trend associated with the decline in flux of AE2g. These two consecutive events have 
an amplitude of about 100 mmag and a duration of 50-150 d, and similar fluctuations BElg and BE2g are evident in the B record. 

Firstly, we analyse the twin events AElg-BElg and AE2g-BE2g. The S /N values for them (the ratios between their semi- 
amplitudes and their mean photometric eiTors) are (S /N)AE\g ~ 4 and {S /N)BEig ~ (S /N)BE2g ~ 3.4. In spite of the fact that AE2g 
is a prominent event, it is poorly traced as a consequence of the 50-day gap and the noisy right wing. Thus, we are not able to 
determine a reliable value of (S /N)AE2g, and the effective signal-to-noise ratio for AE2g could be significantly less than 3-4. The 
difficulties in inferring a time delay from the AE2g-BE2g twin events confirm our suspicions. Unfortunately, it is not possible to 
measure two independent delays, one from AElg-BElg and the other from AE2g-BE2g. The only options are the estimation of a 
delay related to the two flares in the source of intrinsic variability, i.e., using all events in Fig.|6] or a delay corresponding to the first 
flare, i.e., from AElg-BElg. 

Secondly, to calculate the two-flare time delay and magnitude offset (i. e., a constant magnitude shift between the light curves of 
the two quasar compone nts), we use two te chniques: x~ minimization (e.g.. lKundic et al. 119971; lUUan et al. l2006l) and the minimum 
dispersion (D^) method jPelt et al. 111994 [l996.) . characterized by a bin semisize (a) and a decorrelation length (6). The choice of 
Q' = <5 = 9disa good compromise between the A-B connection and time resolution. Through the minimization {a - 9 d), we 
obtain the best solutions of the delay and magnitude offset: AIba = 417 d and Ahiba = -0.090 mag (x^ ~ 1.2). The sign "-" in the 
AniBA value means that the A component is fainter. The minimization (5 = 9 d) gives the best solutions of AIba = 416 d and 
AniBA = -0.092 mag. The uncertainties in the magnitude offset and time delay are inferred from 1000 repetitions of the experiment 
(synthetic light curves based on the observed records). The Icr intervals appear in Table[T] Table [1] indicates that the eiTor in time 
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Fig. 7. Overlapping periods and difference light curves in the g band. We show the overlap between the A (filled circles) and B 
(open circles) whole records, when the A magnitudes are shifted by the best solutions of the time delay and the magnitude offset 
(left panels). We also draw the difference light curve (right panels). The three overlap periods cover ~ 20 d (top panels), ~ 90 d 
(middle panels), and ~ 60 d (bottom panels). 



delay from the;t'^ minimization is substantially less than the error from the minimum dispersion method. Both measurements of the 
two-flare time delay are consistent with the APO main delay in the g band (see Introduction). 

Thirdly, we exclusively use the AElg-BElg twin events. The idea is to accurately measure the gravitational lens delay associated 
with only one flare produced in the source of variability. This time we focus on the 6^ method (see, e.g.. Paper I) and the 
minimization, which produces a delay error smaller than the delay uncertainty from the minimum dispersion technique (see Table[TJ. 
The 6^ technique obtains th e optimal match between the time-shifted discrete autocorrelation function (DAF) and the discrete cross- 
correlation function (DCF: lEdelson & Kroliklll988l) . From the;i'^ minimization (a -9 d), the best solutions of the time delay and 
magnitude offset are 417 d and -0.083 mag, respectively (x^ ~ 0.9). From the 6^ method and 1000 synthetic light curves, the delay 
measurement (Icr interval) is identical to that derived from the technique (see TablelTJ. Therefore, the LT first t win events are 
useful to determine a robust time delay AIba - 417 + 2 d (Icr). This is fully consistent with the APO main delay (Kundi c et al. I 
11997; Goicoechea .2002) . The 6^ analysis also indicates that AIba ^ 424 d (99% confidence interval), so the AElg-BElg delay is 
inconsistent (at about the 3cr level) with the APO secondary delay jGoicoechea Il20 02). The r-band curves of AEl-BEl are not used 
to determine a time delay because S /N <3 for these twin events (e.g., i Piipers Ill997i) . 

We now check for the possible existence of extrinsic variability in our records. We compute the difference light curve between the 
A and B components, since no extrinsic variability should result in a flat difference light curve. To obtain the difference light curve, 
the magnitude- and time-shift ed light curve of image A is subtracted from the light curve of image B (e.g.. lSchmidt & Wambsganssl 
ll998tlGU -Merino et al. Il200"l]) . In Fig.|7](left panels), we show the overlap between the A (filled circles) and B (open circles) whole 
records, when the A magnitudes are shifted by the best solutions of the time delay and the magnitude offset. The difference light 
curve is also plotted in the right panels of Fig.|7] The overlap between A-first season and B-second season covers a very short period 
of about 20 d (see the top panels of Fig.|7]). Fo r this overlap period, the difference curve contains two consecutive deviations from 
the zero level, which are not significative (e.g.. lGil-Merino et al. Il200l1) . The overlap between A-second season and B-third season 
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Fig. 8. Comparison between the DCF (filled circles) and the (DAF) (open circles). While the DCF is the gr cross-correlation 
function, the (DAF) is the average of the gg and rr autocorrelation functions. We use the AE3§-AE3r events (see main text) and 
three bin semisizes: a = 20 (top panel), 25 (middle panel), and 30 (bottom panel) d. 



is much more important than the first overlap (in the top panels). In the middle panels of Fig.|7] we display the situation before the 
50-day gap (see above and Fig.|6]), where the difference curve has a noisy trend that is consistent with zero. In the bottom panels, 
the behaviour after the 50-day gap is shown. In this last period, the difference curve is also mainly noise. However, a clear event 
appears at the beginning of the overlapping period, i.e., six consecutive points are placed above the zero level. Although this naively 
could be interpreted as the wing of a microlensing event (i.e., extrinsic variability), the A data were obtained at the end of a hard 
winter in which the colour coefficient strongly deviated (40-70%) from its average value. While some frames with extreme colour 
coefficients are not considered in the analysis (see the triangles and the open circle around day 3770 in Fig.©, additional adjacent 
frames are also unsuitable for fine variability studies. Therefore, bad weather and anomalous behaviour of the LT devices are the 
most reasonable explanations for the anomalous variation in A that is simultaneously observed in both components. In summary, 
we do not find evidence of extrinsic variability in the light curves of Q0957 -1-561. 



3.2. Delay between optical bands 

The time delay between optical-UV continuum flux variations at two different wavelengths can be used to test the variability 
secenario (e.g.. Collier et al. 1999). It might be produced by reprocessing of high energy radiation in an accretion disc around a 
supermassive black hole. The reprocessing hypothesis assumes that the optical-UV variations are the response of the gas in the 
disc to higher-energy fluctuations in the vecinity of the disc axis. Moreover, the existence of a radiative couplin g between the 
variat i ons is also assumed, i.e., the time delay represents a light-travel time between two disc annuli (see details in ICoUier et al. I 
Il999l) . ICoUier et al. I (1 19991) measured two time lags between fluctuations at two optical wavelengths and the corresponding UV 
fluctuations (UV variabiUty leading optical variations) in the records of NGC 7469 at z - 0.016, and the y fo und a good agr e emen t 
between their delay estimates (~ 1-2 d) and reverberation within an accretion disc. lSergeev et al."] ( 2005l) andlCackett et al. I(l2007b 
also explored the thermal reprocessing hypothesis in local active galactic nuclei. For O0957-H56 1 . IColher reported a delay 
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Fig. 9. Normalised 6^ function from the AE3g-AE3r events. We use bin semisizes a = 20 (dotted line), 25 (dashed line), and 30 
(solid line) d. In Fig. [8j we can observe the presence of time shifts between the DCF and (DAF), which translate into interband 
delay peaks centered on 3-4.5 d (AE3g leading AE3r). 

Table 2. Time lag measurements from the AE3g-AE3r events. 





Time lag^ (d) 


Probability of lags < (%) 


20 


3.0 ± 2.0 


11.6 


25 


4.5 ± 2.5 


6.9 


30 


4.0 ± 2.0 


7.5 


35 


3.5 ± 2.0 


8.5 


40 


3.5 ± 2.0 


9.3 



" We use the technique (see main text) and five values of the bin semisize a. 

* All measurements are lir intervals, and positive lags mean that the r-band event is delayed in relation to the arrival of the associated g-band 
event. 



of about 3.4 d between the r-band and g-band APO main events (g-band events leading those in the r band), which translates into 
a rest-frame lag of about 1.4 d, in excellent agreement with predictions of the disc reprocessing scenario. This first delay betwee n 
optical bands for a GLQ requires an independent confirmation as well as new efforts with other GLQs (e.g.. lKoptelova et al. Il2006h . 
and here we try to reach the first goal. 

For such a task, we focus on the LT events with highest S/N. The AEl-BEl twin events are ruled out because (S /N) < 3 
in the r band. However, there are two very prominent variations around day 4150 in the top panels of Figs. |4l-|5](A component). 
These AE3g and AE3r variations last ~ 250 d (the whole light curves of A in the third season are considered as large events) and 
have signal-to-noise ratios above 6. We use fluxes in arbitrary units / = 10^ x lO"'^*" to compare AE3g and AE3r. The use of 
fluxes (instead of magnitudes m) permits a fair cross-coiTelation between two records that, apart from a possible delay, differ in a 
multiplicative constant and an additive constant. On average, the light curves were sampled two times per week. However, there are 
20-day gaps around day 4180. Unfortunately, due to a combination of the kind of variability (time asymmetric events consisting of 
slow rises and rapid declines) and these short gaps close to the maxima, it is not possible to infer a reliable DCF with good time 
resolution, i.e., a < 10 d (see above). For a < 20 d, 20-day artifacts at lags of + 50 d appear in the DCF. This unphysical signal 
at + 50 d is only avoided using longer bins with a > 20 d, so we are forced to take relatively long bins. This is not a problem 
at all, but the measurement would be more accurate (but not more reliable) with better time resolution. Some DCF (filled circles) 
and (DAF) (open circles) trends are shown in Fig. [8] The top, middle, and bottom panels of Fig. [8] contain the results for a = 20, 
25, and 30 d, respectively. Here, {DAF} is the average of the gg and rr autocorrelation functions, whereas DCF represents the gr 
cross-correlation function. 

In Fig. [8] there are no important distortions in the features of the DCF compared to the features in the (DAF), but the existence 
of a delay of several days is evident. In other words, to get an optimal match, the {DAF} should be shifted by several days. Possible 
values of this time shift (0) versus the a ssociated 6^(0) values nor malised by its minimum value 6^(0o) are plotted in Fig. |9] The 
6^(6) function was defined in Eq. (7) of ISerra-Ricart et alTI d 19991) (see also above), and we use a = 20 (dotted line), 25 (dashed 
line), and 30 (solid line) d in Fig.|9l This figure displays relatively narrow peaks centered on 3-4.5 d (best values of the interband 
delay; AE3g leading AE3r). Uncertainties are again computed by applying the 6^ minimization to 1000 synthetic data sets. Through 
the distributions of delays (a = 20-40 d), five Icr measurements are presented in Table |2l The 6^ results in Table |2] agree with the 
previous time lag determination from APO light curves, and we adopt Af^g - 4.0 + 2.0 d (using an intermediate bin semisize a = 30 
d; the probability of Atrg < is only 7.5%) as our final Icr measurement. 
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4. Summary and conclusions 

Liverpool Quasar Lens Mo nitoring is a long-t erm project to follow the optical (gri bands) variability of 10-20 GLQs with the 
Liverpool robotic telescope (ISteele et al. II2004I) . The first phase of this project (LQLM I) was conducted between January 2005 and 
July 2007. While in Paper I we mainly studied the intrinsic variability of Q0909+532 in the r band, in this paper we present the 
monitoring programme of Q0957-I-561 in the gr bands. A main goal of our project (LQLM) is to considerably increase the public 
database of GLQs. Thus, all LQLM I pre-processed frames of Q0909+532 and Q0957-I-561 are publicly available on the Lens Image 
Archive of the German Astrophysical Virtual Observatorjfl 

We have fully developed two photometric pipelines through the 3 years of observations and analyses. The transformation 
pipeline incorporates zero-point, colour, and inhomogeneity corrections in the instrumental fluxes, so photometry to the 1-2% 
is achieved for Q0957-I-561A and Q0957+561B. We detect an inhomogeneous response over the CCD area, whi ch has an amph- 
tude of ~ 80 mmag (from maximum to minimum) and is consistent with studies in other optical telescopes (e.g.. iManfroid et aTl 
120011: iMagnier & Cuillandre 1 12004 . Moreover, the colour coeflicient is allowed to vary through time, because the atmospheric- 
instrumental conditions signicantly evolve through 2.5 years of monitoring. Due to atmospheric-instrumental problems at some 
epochs, the colour coefficient reaches anomalous values, i.e., we obtain dramatic deviations with respect to the average coefficient. 
Thus, the frames corresponding to an anomalous coefficient are removed or not considered. 

The LT gr light curves of Q0957+561 show several prominent events and gradients, and some of them (in the g band) are used to 
infer a time delay between components AIba = 417 + 2 d (Icr). This gravitationa l lens delay from new ^-band events is in agreement 
with the delay from the previous APO g-band main events (Kundic et al. 1997), so the associated UV flares in the variable source 
(APO and LT events) probably originate in the same emission region ( Y onehara 1999) . Ta king into accou nt that the previous APO 
^r-band main events are plausibly due to reverberation within an irradiated accretion disc (ICoUier 11200 lb . the new gr-hand events 
are likely related to flares in the central accretion disc. In addition, the delay between the t wo new LT larg e events in the g and r 
bands: Atrg - 4.0 + 2.0 d (Icr; the ^-band event is leading), coincides with the estimation bv lCoUierl (12001 1) and agrees with flares 
generated during reprocessing in the accretion disc. Therefore, most APO-LT variations in the g and r bands are very probably 
associated with the gas disc around the supermassive black hole (only the APO secondary events have been associated with flares 
that were produced far away from the accretion disc; see Introduction). The detection of the same interband delay (between the g 
and r band) in the two monitoring campaigns (APO and LT) also suggests that the accretion disc reprocessing in Q0957-H561 is a 
usual occurance at different times for different prominent flares. Hence, very likely, most observed variations in the g and r bands 
(APO and LT fluctuations with an amplitude of ~ 100 mmag and lasting ~ 100 d) are associated with reverberation within the gas 
disc around the supermassive black hole. 

We add 2.5 years of time coverage to the previous 1. 5-year gr-band records of Q0957-I-561, and remark that our difference light 
curves are consistent with zero. Thus, there is no evidence of extrinsic variations in both APO and LT independent experiments 
separated by ~ 10 years. These results disagree with the claim by Schild (2005) that microlenses in the lensing galaxy affect 
the observed variability. Therefore, the complex quasar structure suggested by this author is not supported by the gr-hand light 
curves of Q0957-I-561. We also remark that double replicas in the records of the B component are not detected in the APO and LT 
experiments. This clearly contradicts previous conclusions bv lHirv et al. I ( l2007h . which indicated that the full B light curve of the 
lensed quasar can be decomposed into a sum of two similar and time shifted curves. Finally, the APO-LT combined database of 
Q0957 -1-561 (together with other monitorings done between both experiments) is a promising t ool for studying the quasar structure 
and the composition of the lensing halo (e.g.. lSchmidt & Wambsganss lll998l ; lKochanekll2004l) . 
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Appendix A: Transformation Equations 

The initial transformation equations for a given reference star are 

g*{tj) = g + Agitj) + Cg{tj){g - r), (A.l) 
r\tu) = r + A,{t,) + C,.{tu){r-i), (A.2) 

where g* and r* are the instrumental magnitudes of the star, g, r, and / are its standard magnitudes, Ag and A^ are the zero-point 
terms (including instrumental and atmospheric corrections), and Cg and Cr are the colour coefficients. The zero-point terms and 
the colour coefficients are expected to significantly change during the 2.5-year monitoring period, so we explicitly consider their 
time evolution. Here, tj and tk denote observation times in the g and r bands, respectively. As we are initially interested in the usual 
systematic corrections, Eqs. (A. 1-2) do not include other possible terms (see here below). Instead of the LT photometric system 
(ugriz = u'g' r'i'z'), we want to use the SPSS "natural" system, since accurate standard magnitudes are available in this SDSS 2.5m 
system (e.g.. Smi th et al. II2002I : IStoughto n et al. II2002I) . SDSS magnitudes are also suitable for comparing our results with future 
data of Q0957-H561 using different facilities and/or SDSS quasar studies/databases (e.g., Vand en Berk et al. .2004; Schneider et al. I 
|2007|) . From equations for transforming LT magnitudes to magnitudes in the SDSS systerrQ: 

g = gsDss + Bgig -r) + Kg, (A.3) 
= rsDss + Brir - i) + K,-, (A.4) 

and equations that relate LT and SDSS colours: 

g-r ^ Ogrig - r)sDss + bgr, (A.5) 
r - / = a,i(r - i)sDss + byi, (A.6) 

it is possible to rewrite Eqs. (A. 1-2) as 

g\tj) = gsDss + agitj) +/3gitj)ig - r)sDss, (A.7) 

r'itk) = rsDss +ar(tk)+Mtk)(r-i)sDss- (A.8) 
The Og term and the /3g coefficient are given by (it is trivial to write expressions for a,- and /3r) 

agitj) = Agitj) + Kg + bgABg + Cgitj)], (A.9) 

Pgitj) = agABg + Cgitj)]. (A.IO) 



^ See the Web site |http ://www. sdss . org/dr6/algorithms71 
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Table A.l. Adopted standard magnitudes of the reference stars. 



Star 


gSDSS 


I'SDSS 


isDSS 


X 


14.213 


13.849 


13.750 


G 


14.461 


14.157 


14.060 


F 


14.513 


14.186 


14.089 


H 


15.116 


14.422 


14.174 


D 


15.485 


14.95 P 


14.770 


E 


15.816 


15.21~ 


15.018 


R 


17.879 


16.801 


16.419 



" The SDSS catalogue seems to contain a wrong value of the r-SDSS magnitude of the D star (rsoss = 15.674), so the D star would be fainter 
than the E star in this band. This disagrees with our current LT observations and several previous observations in the red region of the optical 
spectrum. Thus, the r-SDSS magnitude is inferred through the rguss vs. VR relationship: rgoss = V - 0.89(K - R) + 0.3 9. This law is based on 
the r-SDSS magnitudes of the rest of stars and the corresponding VR magnitudes in Tables 1-2 of lOvaldsen et al. I ( l2003bh . 



Taking into account typical values of Q ( 0.029) and C,- (~ 0.034) reported by the LT team (on the LT Web site; see main text), 

the SDSS estimates of Bg (~ -0.060) and Br (~ -0.035) on the SDSS Web site (see here above), and Ogr ~ an ~ 1, we expect 

typical colour coefficients fig 0.089 and yS,- ~ -0.001. On the other hand, the adopted standard magnitudes of the reference stars 

(XGFHDER field stars; see main text) appear in Table I A. ll These are PSF magnitudes in the SDSS catalogue^. 

In order to achieve 1-2% photometric accuracy with the RATCam camera (on the LT), one additional detail must be taken 
into account in the transformation equations (A. 7-8). We introduce an inhomogen eity term that corrects the flat-field systematic 
error over the camera area, which might have a total amplitude of ~ 50 mmag (e.g.. iManfroid et ani200ll; iMagnier & Cuillandre I 
|2004|) . For example, this kind of error could be related to twilight flats. During twilight exposures, some scattered light (within the 
dome) would reach the camera, and thus, the illumination would not be homogeneous. This effect invalidates the basic hypothesis 
of homogeneous illumination. Here, we assume a second order 2D polynomial to account for the inhomogeneity term, so the final 
transformation equations are 

g'itj) ^ gSDSS +ag{tj)+Pg{tj)(g-r)sDSS + ^ yg,nm^{tj)y"'(tj), (A. 11) 

0<«+m<2 

r\tk) ^ rsDSS +ar{tk)+l3r{tk){r-i)sDSS + ^ yr,nmX"{tk)y"'{tk), (A.12) 

0<«+m<2 

where {x,y) is the 2D position of the star on the CCD. To find the relevant parameters in the g band, i.e., agitj), /Sgitj), and jg.nm, 
we may fit the observed magnitude deviations (instrumental - standard) of the seven reference stars to the model incorporating the 
three systematic terms: zero-point, colour, and inhomogeneity. Once the fit has been made, the g-SDSS magnitude of any point-like 
source (star or quasar) is derived in a straightforward way: 

giSDSS) = gsDss g\tj) - agitj) -I3g{tj){g - r)sDss - J] 7g.nmx"{tj)f'{tj). (A.13) 

0<n+m<2 

In Eq. (A.13), 5 represents the deviation caused by random noise (e.g., photon noise) and unkown (but presumibly small) systematic 
corrections. For a non-variable star (e.g., a reference star), variations in g{S DSS) are generated by noise (S). However, for variable 
stars or quasars, there are two kinds of variability. While true variability is due to time evolution of gsoss, noise is an additional 
source of fluctuations. The (g - r)sDss colour of Q0957-I-561 might also evolve over time. Thus, the use of an average colour 
(ig - t')sDss) in the colour correction introduces a systematic noise Smi = PMd^ig ~ Qspss associated with the colour variation. 
Fortunately, for moderate fluctuations with amplitude of ~ 25 mmag (e.g.. lKundic et al.1l 19951) . the amplitude of the colour noise is 
only ~ 0.2%, well below our accuracy goal. The r-SDSS magnitude of a source is given by an expression similar to Eq. (A. 13). 
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